%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cs\_mass\_flux}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Le but de ce sous-programme est principalement de calculer le flux de masse aux
faces. Il prend une variable vectorielle associée au centre des cellules
(généralement la vitesse), la projette aux faces en la multipliant par la
masse volumique, et la multiplie scalairement par le vecteur surface.
Plus généralement, \fort{cs\_mass\_flux} est aussi appelé comme première étape
du calcul d'une divergence (terme en $\dive(\rho\tens{R})$ en
$R_{ij}-\varepsilon$, filtre Rhie \& Chow, ...).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discrétisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

La figure \ref{Base_Inimas_fig_geom} rappelle les diverses définitions géométriques
pour les faces internes et les faces de bord. On notera
$\displaystyle \alpha=\frac{\overline{FJ^\prime}}{\overline{I^\prime J^\prime}}$ (défini aux faces
internes uniquement).

\begin{figure}[h]
\parbox{8cm}{%
\centerline{\includegraphics[height=4cm]{facette}}}
\parbox{8cm}{%
\centerline{\includegraphics[height=4cm]{facebord}}}
\caption{\label{Base_Inimas_fig_geom}Définition des différentes entités
géométriques pour les faces internes (gauche) et de bord (droite).}
\end{figure}


\subsection*{Faces internes}
On ne conna\^\i t pas la masse volumique à la face, cette dernière doit donc
aussi être interpolée. On utilise la discrétisation suivante :

\begin{equation}
(\rho \vect{u})_F = \alpha(\rho_I \vect{u}_I)
+(1-\alpha)(\rho_J \vect{u}_J)
+\ggrad\!(\rho\vect{u})_O.\vect{OF}
\end{equation}
La partie en $\alpha(\rho_I \vect{u}_I)
+(1-\alpha)(\rho_J \vect{u}_J)$ correspondant en fait à
$(\rho\vect{u})_O$. Le gradient en $O$ est calculé par interpolation :
$\displaystyle\ggrad\!(\rho\vect{u})_O=
\frac{1}{2}\left[\ggrad\!(\rho\vect{u})_I+\ggrad\!(\rho\vect{u})_J\right]$. La
valeur $\displaystyle\frac{1}{2}$ s'est imposée de manière heuristique au
fil des tests
comme apportant plus de stabilité à l'algorithme global qu'une interpolation
faisant intervenir $\alpha$. L'erreur commise sur $\rho\vect{u}$ est en
$O(h^2)$.


\subsection*{Faces de bord}
Le traitement des faces de bord est nécessaire pour y calculer le flux de
masse, bien sûr, mais aussi pour obtenir des conditions aux limites pour le
calcul du $\ggrad\!(\rho \vect{u})$ utilisé pour les faces internes.

Pour les faces de bord, on conna\^\i t la valeur de $\rho_F$, qui est stockée
dans la variable \var{ROMB}. De plus, les conditions aux limites pour $\vect{u}$
sont données par des coefficients $A$ et $B$ tels que :
\begin{equation}
u_{k,F} = A_k + B_ku_{k,I^\prime} =
A_k + B_k\left(u_{k,I} + \grad\!(u_k)_I.\vect{II^\prime}\right)
\end{equation}
($k\in\{1,2,3\}$ est la composante de la vitesse, l'erreur est en $O(B_kh)$)

On a donc à l'ordre 1 :
\begin{equation}
(\rho u_k)_F = \rho_F\left[A_k + B_k\left(u_{k,I} +
\grad\!(u_k)_I.\vect{II^\prime}\right)\right]
\end{equation}

Mais pour utiliser cette formule, il faudrait calculer $\ggrad\!(\vect{u})$ (trois
appels à \fort{GRDCEL}), alors qu'on a déjà calculé
$\ggrad\!(\rho\vect{u})$ pour les faces internes. Le surcoût en temps serait alors
important. On réécrit donc :
\begin{eqnarray}
(\rho u_k)_F & = & \rho_F A_k + \rho_F B_ku_{k,I^\prime}\\
& = & \rho_F A_k + B_k\frac{\rho_F}{\rho_{I^\prime}}(\rho u_k)_{I^\prime}
\label{Base_Inimas_eq_rhoufacea}\\
& = & \rho_F A_k + B_k\frac{\rho_F}{\rho_{I^\prime}}(\rho u_k)_I
+B_k\frac{\rho_F}{\rho_{I^\prime}}\grad\!(\rho u_k)_I.\vect{II^\prime}
\label{Base_Inimas_eq_rhoufaceb}
\end{eqnarray}

Pour calculer les gradients de $\rho\vect{u}$, il faudrait donc en théorie
utiliser les coefficients de conditions aux limites équivalents :\\
$\tilde{A}_k = \rho_F A_k$\\
$\displaystyle \tilde{B}_k = B_k\frac{\rho_F}{\rho_{I^\prime}}$

Ceci paraît délicat, à cause du terme en
$\displaystyle \frac{\rho_F}{\rho_{I^\prime}}$, et en particulier à l'erreur
que l'on peut commettre sur $\rho_{I^\prime}$ si la reconstruction des gradients
est imparfaite (sur des maillages fortement non orthogonaux par exemple).
On réécrit donc l'équation
(\ref{Base_Inimas_eq_rhoufaceb}) sous la forme suivante :
\begin{equation}
(\rho u_k)_F=\rho_F A_k + B_k\frac{\rho_I\rho_F}{\rho_{I^\prime}}u_{k,I}
+B_k\frac{\rho_F}{\rho_{I^\prime}}\grad\!(\rho u_k)_I.\vect{II^\prime}
\end{equation}


Pour le calcul du flux de masse au bord, on va faire deux approximations. Pour
le deuxième terme, on va supposer $\rho_{I^\prime}\approx\rho_I$ (ce qui
conduit à une erreur en $O(B_kh)$ sur $\rho\vect{u}$ si
$\rho_{I^\prime}\ne \rho_I$). Pour le
troisième terme, on va supposer $\rho_{I^\prime}\approx\rho_F$. Cette
dernière approximation est plus forte, mais elle n'intervient que dans la
reconstruction des non-orthogonalités ; l'erreur finale reste donc faible
(erreur en $O(B_kh^2)$ sur $\rho\vect{u}$ si
$\rho_{I^\prime}\ne \rho_F$).
Et au final, le flux de masse au bord est calculé par :
\begin{equation}
\dot{m}_F = \sum\limits_{k=1}^{3}\left[\rho_F A_k + B_k\rho_Fu_{k,I}
+B_k\grad\!(\rho u_k)_I.\vect{II^\prime}\right]S_k
\end{equation}

Pour le calcul des gradients, on repart de l'équation (\ref{Base_Inimas_eq_rhoufacea}), sur
laquelle on fait l'hypothèse que $\rho_{I^\prime}\approx\rho_F$. Encore une
fois, cette hypothèse peut être assez forte, mais les gradients obtenus ne
sont utilisés que pour des reconstructions de non-orthogonalités ; l'erreur
finale reste donc là encore assez faible.
Au final, les gradients sont calculés à partir de la formule suivante :
\begin{equation}
(\rho u_k)_F = \rho_F A_k + B_k(\rho u_k)_{I^\prime}
\end{equation}
ce qui revient à utiliser les conditions aux limites suivantes pour
$\rho \vect{u}$:\\
$\tilde{A}_k = \rho_F A_k$\\
$\tilde{B}_k = B_k$

\minititre{Remarque}

Dans la plupart des cas, les approximations effectuées n'engendrent aucune
erreur. En effet :\\
- dans le cas d'une entrée on a généralement $B_k=0$, avec un flux de
masse imposé par la condition à la limite.\\
- dans le cas d'une sortie, on a généralement flux nul sur les scalaires
donc sur $\rho$, soit \mbox{$\rho_F=\rho_{I^\prime}=\rho_I$}.\\
- dans le cas d'une paroi, on a généralement $B_k=0$ et le flux de masse
est imposé nul.\\
- dans le cas d'une symétrie, on a généralement
$\rho_F=\rho_{I^\prime}=\rho_I$ et le flux de masse est imposé nul.\\
Pour sentir un effet de ces approximations, il faudrait par exemple une paroi
glissante ($B_k\ne0$) avec un gradient de température ($\rho_F\ne\rho_I$).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

La vitesse est passée par les arguments \var{UX}, \var{UY} et \var{UZ}. Les
conditions aux limites de la vitesse sont \var{COEFAX}, \var{COEFBX}, ... Le
flux de masse résultat est stocké dans les variables \var{FLUMAS} (faces
internes) et \var{FLUMAB} (faces de bord). \var{QDMX}, \var{QDMY} et \var{QDMZ}
sont des variables de travail qui serviront à stocker $\rho\vect{u}$, et
\var{COEFQA} servira à stocker les $\tilde{A}$.

\etape{Initialisation éventuelle du flux de masse}
Si \var{INIT} vaut 1, le flux de masse est remis à zéro. Sinon, le
sous-programme rajoute aux variables \var{FLUMAS} et \var{FLUMAB} existantes le
flux de masse calculé.


\etape{Remplissage des tableaux de travail}
$\rho\vect{u}$ est stocké dans \var{QDM}, et $\tilde{A}$ dans \var{COEFQA}.


\etape{Cas sans reconstruction}
On calcule alors directement\\
$\displaystyle \var{FLUMAS}=\sum\limits_{k=1}^{3}\left[
\alpha(\rho_I u_{k,I})+(1-\alpha)(\rho_J u_{k,J})\right]S_k$\\
et\\
$\displaystyle \var{FLUMAB}=\sum\limits_{k=1}^{3}\left[
\rho_F A_k + B_k\rho_Fu_{k,I}\right]S_k$


\etape{Cas avec reconstruction}
On répète trois fois de suite les opérations suivantes, pour $k=1$, 2 et 3
:\\
- Appel de \fort{GRDCEL} pour le calcul de $\grad\!(\rho u_k)$.\\
- Mise à jour du flux de masse\\
$\displaystyle \var{FLUMAS}=\var{FLUMAS} + \left[
\alpha(\rho_I u_{k,I})+(1-\alpha)(\rho_J u_{k,J})
+\frac{1}{2}\left[\grad\!(\rho u_k)_I+\grad\!(\rho u_k)_J\right]
.\vect{OF}\right]S_k$\\
et\\
$\displaystyle \var{FLUMAB}=\var{FLUMAB}+\left[
\rho_F A_k + B_k\rho_Fu_{k,I}
+B_k\grad\!(\rho u_k)_I.\vect{II^\prime}\right]S_k$


\etape{Annulation du flux de masse au bord}
Quand le sous-programme a été appelé avec la valeur \var{IFLMB0=1}
(c'est-à-dire quand il est réellement appelé pour calculer un flux de
masse, et pas pour calculer le terme en $\dive(\rho\tens{R})$ par exemple), le flux
de masse au bord \var{FLUMAB} est forcé à 0, pour les faces de paroi et de
symétrie (identifiées par \var{ISYMPA=0}).
